#include "matrix.h"

Matrix::~Matrix()
{
    delete []datap;
}

Matrix::Matrix(const Matrix &a)
{
    row = a.row;
    col = a.col;
    datap = (double *)malloc(row * col * sizeof(double));
    for (int i = 0; i < row * col; i++)
    {
        datap[i] = a.datap[i];
    }
}

Matrix::Matrix(int _row, int _col, int cvalue)
{
    row = _row;
    col = _col;
    datap = (double *)malloc(row*col*sizeof(double));
    for (int i = 0; i < row * col; i++)
    {
        datap[i] = cvalue;
    }
}

Matrix::Matrix(int _row, int _col, double *p)
{
    row = _row;
    col = _col;
    datap = (double *)malloc(row*col*sizeof(double));
    for (int i = 0; i < row * col; i++)
    {
        datap[i] = p[i];
    }
}

void Matrix::show()
{
    for (int i = 0; i < row; i++)
    {
        for (int j = 0; j < col; j++)
        {
            printf("%f ", datap[i*col + j]);
        }
        printf("\n");
    }
    printf("\n");
}

Matrix Matrix::diag()
{
    Matrix a(row, 1, 0);
    for (int i = 0; i < row; i++)
    {
        a.datap[i] = datap[i * (col + 1)];
    }
    return a;
}

void Matrix::trans()
{
    double *q = (double *)malloc(row*col*sizeof(double));
    for (int j = 0; j < col; j++)
    {
        for (int i = 0; i < row; i++)
        {
            q[i + j * row ] = datap[i * col + j];
        }
    }
    for (int k = 0; k < row * col; k++)
    {
        datap[k] = q[k];
    }
    std::swap(row, col);
}

Matrix Matrix::cut(int rowbegin, int rowend, int colbegin, int colend)
{
    Matrix a(rowend - rowbegin + 1, colend - colbegin + 1, 0);
    for (int i = rowbegin; i < rowend + 1; i++)
    {
        for (int j = colbegin; j < colend + 1; j++)
        {
            a.datap[(i - rowbegin) * a.col + (j - colbegin)] = datap[(i - 1) * col + j - 1];
        }
    }
    return a;
}

double Matrix::norm0()
{
    double max = 0;
    double sum;
    for (int i = 0; i < row; i++)
    {
        sum = 0;
        for (int j = 0; j < col; j++)
        {
            sum += fabs(datap[i * col + j]);
        }
        if (sum>max)
        {
            max = sum;
        }
    }
    return max;
}

double Matrix::norm1()
{
    double max = 0;
    double sum;
    for (int i = 0; i < col; i++)
    {
        sum = 0;
        for (int j = 0; j < row; j++)
        {
            sum += fabs(datap[j * col + i]);
        }
        if (sum>max)
        {
            max = sum;
        }
    }
    return max;
}

double Matrix::norm()
{
    Matrix a = (*this);
    Matrix b = a;
    a.trans();
    b = a * b;
    Matrix u(b.row, 1, 1);
    double m = -1;
    double eps = 1e-5;
    Matrix y(b.row, 1, 1);
    while (vnorm(b*u-u*m)>eps)
    {
        y = b * u;
        m = y.norm0();
        u = y * (1 / m);
    }
    return sqrt(m);
}

Matrix Matrix::idown()
{
    Matrix a = ones(row, col);
    for (int i = 0; i < row; i++)
    {
        for (int j = 0; j <= i; j++)
        {
            if (i!=j)
                a.datap[i*a.col+j] = datap[i*col+j];
            else
                a.datap[i * a.col + j] = 1;
        }
    }
    return a;
}

Matrix Matrix::down()
{
    Matrix a = ones(row);
    for (int i = 0; i < row; i++)
    {
        for (int j = 0; j <= i; j++)
        {
            a.datap[i*a.col+j] = datap[i*col+j];
        }
    }
    return a;
}

Matrix Matrix::downarc()
{
    Matrix a = ones(row);
    for (int j = 1; j < col; j++)
    {
        for (int i = j+1; i < row + 1; i++)
        {
            a.datap[(i - 1) * col + j - 1] = -datap[(i - 1) * col + j - 1]/datap[(j - 1) * (col + 1)];
            for (int m = 1; m < j; m++)
            {
                a.datap[(i - 1) * col + m - 1] += a.datap[(j - 1) * col + m - 1] * a.datap[(i - 1) * col + j - 1];
            }
        }
        for (int k = j+1; k < col + 1; k++)
        {
            for (int l = j+1; l < row+1; l++)
            {
                datap[(l - 1) * col + k - 1] -= datap[(l - 1) * col + j - 1] * datap[(j - 1) * col + k - 1];
            }
        }
    }
    for (int i = 0; i < row; i++)
    {
        for (int j = 0; j <= i; j++)
        {
            a.datap[i * a.col + j] /= datap[i * (a.col + 1)];
        }
    }
    return a;
}

Matrix Matrix::arc()
{
    Matrix a = ones(row);
    Matrix b = *this;
    double value;
    for (int j = 1; j < b.col; j++)
    {
        for (int i = j+1; i < b.row + 1; i++)
        {
            value = b.datap[(i - 1) * col + j - 1]/b.datap[(j - 1) * (col + 1)];
            for (int m = 1; m < j+1; m++)
            {
                a.datap[(i - 1) * col + m - 1] -= a.datap[(j - 1) * col + m - 1] * value;
            }
            for (int m = 1; m < col+1; m++)
            {
                b.datap[(i - 1) * col + m - 1] -= b.datap[(j - 1) * col + m - 1] * value;
            }
        }   
    }
    for (int i = 0; i < row; i++)
    {
        for (int j = 0; j <= i; j++)
        {
            a.datap[i * a.col + j] /= b.datap[i * (a.col + 1)];
        }
        for (int j = col-1; j >=i; j--)
        {
            b.datap[i * a.col + j] /= b.datap[i * (a.col + 1)];
        }
    }
    for (int j = col; j >1; j--)
    {
        for (int i = j-1; i > 0; i--)
        {
            for (int m = 1; m < col+1; m++)
            {
                a.datap[(i - 1) * col + m - 1] -= b.datap[(i - 1) * col + j - 1] * a.datap[(j - 1) * col + m - 1];
            }
            b.datap[(i - 1) * col + j - 1] -= b.datap[(i - 1) * col + j - 1] * b.datap[(j - 1) * (col +1)];
        }
    }
    return a;
}

Matrix operator*(const Matrix &a, double b)
{
    Matrix c(a.row, a.col, 0);
    for (int i = 0; i < a.row*a.col; i++)
    {
        c.datap[i] = a.datap[i] * b;
    }
    return c;
}

Matrix operator*(const Matrix &a, const Matrix &b)
{
    Matrix c(a.row, b.col, 0);
    for (int i = 0; i < a.row; i++)
    {
        for (int j = 0; j < b.col; j++)
        {
            for (int k = 0; k < a.col; k++)
            {
                c.datap[i * b.col + j] += a.datap[i*a.col+k] * b.datap[j+k*b.col];
            }
        }
    }
    return c;
}

Matrix Matrix::operator=(const Matrix &a)
{
    row = a.row;
    col = a.col;
    for (int i = 0; i < row * col; i++)
    {
        datap[i] = a.datap[i];
    }
    return *this;
}

Matrix Matrix::operator+(const Matrix &a)
{
    double p[a.row * a.col];
    for (int i = 0; i < a.row * a.col; i++)
    {
        p[i] = a.datap[i]+this->datap[i];
    }
    return Matrix(a.row, a.col, p);
}

Matrix Matrix::operator-(const Matrix &a)
{
    double p[a.row * a.col];
    for (int i = 0; i < a.row * a.col; i++)
    {
        p[i] = this->datap[i] - a.datap[i];
    }
    return Matrix(a.row, a.col, p);
}

double vnorm(Matrix x)
{
    double a = 0;
    for (int i = 0; i < x.row; i++)
    {
        a += x.datap[i] * x.datap[i];
    }
    a = sqrt(a);
    return a;
}

Matrix ones(int n)
{
    Matrix a(n, n, 0);
    for (int i = 0; i < n; i++)
    {
        a.datap[i * (a.col + 1)] = 1;
    }
    return a;
}

Matrix ones(int m, int n)
{
    Matrix a(m, n, 0);
    for (int i = 0; i < n && i<m; i++)
    {
        a.datap[i * (a.col + 1)] = 1;
    }
    return a;
}

Matrix diag(int n, double *p)
{
    Matrix a(n, n, 0);
    for (int i = 0; i < n; i++)
    {
        a.datap[i * (a.col + 1)] = p[i];
    }
    return a;
}

Matrix trans(Matrix a)
{
    Matrix b(a);
    b.trans();
    return b;
}

/*int main()
{
    double p[9] = {1, 2, 3, 4, 5, 6, 7, 8, 9};
    Matrix a(3, 3, p);
    a.show();
    Matrix b = a;
    std::cout << std::endl;
    b.trans();
    b.show();
    Matrix c = a * b;
    std::cout << std::endl;
    c.show();
    Matrix d = a.down();
    Matrix c = d.downarc();
    d.show();
    double p[9] = {1, 2, 7, 3, 4, 6, 5, 9, 8};
    Matrix a(3, 3, p);
    Matrix d = a.arc();
    Matrix e = d * a;
    e.show();
    a.show();
    return 0;
}*/